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A MAP OF A PROPERTY 
FIELD OF THE INVENTION 

The present invention relates to a method of creating a 
map of a property of an object and to a map when created 
in accordance with the method. The invention relates 
particularly, though not exclusively, to a method of 
creating a map of a property derived from nuclear 
magnetic resonances related to a line or region within an 
object . 

BACKGROUND OF THE INVENTION 

It is currently of significant interest in different 
areas of research and analysis, in particular in 
medicine, to obtain non-destructively images or line- 
scans related to regions within objects. Typical objects 
may include biological specimen such as live human 
brains. The analysis of such images or line -scans may 
reveal information that is essential for understanding 
the functional organisation of the human brain. 

In order to gain a more detailed understanding of the 
activation of the cerebral cortex of the brain, the 
anatomical organisation of the brain has to be known. 
With regard to the cerebral cortex, which is of 
particular interest in studies of cognition, a key 
question is how one relates the results of a functional 
neuro-imaging study with particular cortical areas in 
which changes in neural activity occur. 

An enormous and costly effort is expended in carrying out 
functional neuro-imaging of ever increasing resolution, 
but in the end the researcher is left with only a crude 
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anatomical estimate of the activated cortical areas. If 
this situation is to be resolved, better human brain 
atlases are required. However, these suffer from one 
major problem; namely inter-subject variations in brain 
micro- and macrostructure . 

Histological examination of the cyto-architecture , myelo- 
architecture and chemo-architecture of the brain remains 
the best method of accurately defining anatomically 
distinct regions. Ideally, if this can be performed in 
every subject used in functional neuroimaging studies, 
each individual's functional results can be precisely and 
accurately related with their own anatomical brain map. 
Furthermore, if a functionally active cortical focus were 
consistently found to be correlated with a well defined 
anatomical region across a population of subjects, this 
would significantly advance understanding of human brain 
organisation. Obviously, this is not practical except in 
a few rare cases. Instead, the results from postmortem 
brains of "non-sub j ects" in functional studies are used 
which again results in problems of inter-subject 
variability. 

SUMMARY OP THE INVENTION 

According to a first aspect of the invention there is 
provided a method of mapping a property of a three 
dimensional object, said method comprising the steps of 

mapping the property for at least a portion of the 

object, 

providing a line or region which defines intersections 
with part of said mapped portion, and 
displaying the property for the intersections. 



WO 02/098292 PCT/AU02/00731 

3 

According to a second aspect of the invention there is 
provided a method of mapping a property of a three 
dimensional object, said method comprising the steps of: 

mapping the property for a plurality of slices within 

the object, 

providing a line or region which defines intersections 
with part of at least some of the plurality of slices, 
and 

displaying the property for the intersections. 

Preferably the property relates to nuclear magnetic 
resonances. More preferably, mapping of the property for 
the slices results in a plurality of nuclear magnetic 
resonance images . 

It should be understood that nuclear magnetic resonances 
and nuclear magnetic resonance images are terms usually 
used by physicists and chemists. It is common practice to 
abbreviate these terms in a clinical environment to 
magnetic resonances and magnetic resonance images, 
respectively . 

Preferably the line or region is a line. More preferably, 
the line is a plurality of lines. 

Preferably, each of the lines is substantially 
perpendicular to the surface of the three-dimensional 
object . 

In one embodiment the object is a biological object such 
as the brain of a subject (or another organ of a human or 
animal, or a semi -biological or inanimate object) and 
each of the lines is a line through the cortex of the 
brain. The cortex of the human brain is convoluted and 
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contains a plurality of grooves. Therefore, a line 
substantially perpendicular to the surface, for example 
within a groove, would typically intersect not only one 
slice, but can intersect a plurality of slices which are 
related to nuclear magnetic resonance images taken from 
adjacent areas of the brain. 

Preferably, the method further comprises the steps of 
approximating the shape of the object by a three 
dimensional lattice of two-dimensional forms created by a 
first computer routine. More preferably, the lattice is a 
mesh. 

Preferably, the lines are each substantially 
perpendicular to a respective one of the two-dimensional 
forms . 

Preferably the intersections of each of the lines with 
each of the slices are selected using a second computer 
routine. More preferably a third computer routine maps 
nuclear magnetic resonance values for the intersections 
by displaying an array of values . 

Preferably the method further comprises the step of 
analysing the array of values using a fourth computer 
routine. More preferably the fourth computer routine 
analyses quantitative properties. 

Preferably analysing the array of values using the fourth 
computer routine comprises differentiation. More 
preferably analysing the array of values comprises 
determining the positions of the maxima, minima and/or 
zero -points of the array of values and the derivatives of 
the array of values . 
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In one embodiment, the positions of the inner and outer 
boundaries of the grey matter of the brain may be 
approximated by determining the positions of the maxima 
of the first derivative of the array of values. Two 
closely followed zero-points of the first derivative may 
indicate the presence of lamination. By measuring the 
relative positions of the zero points the thickness of 
the lamination can be determined. 

According to a third aspect of the invention there is 
provided a map produced by a method of mapping a property 
of a three dimensional object, said method comprising the 
steps of : 

mapping the property for at least a portion of the 
object , 

providing a line or region which defines intersections 
with part of said mapped portion, and 
displaying the property for the intersections. 

According to a fourth aspect of the invention there is 
provided a map produced by a method of mapping a property 
of a three dimensional object, said method comprising the 
steps of : 

mapping the property for a plurality of slices within 
the object, 

providing a line or region which defines intersections 
with part of at least some of the plurality of slices, 
and 

displaying the property for the intersections. 

According to a fifth aspect of the invention there is 
provided a method of mapping a property of an object, 
said method comprising the steps of : 
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mapping the property for a sxice within the object, 
providing a line within the slice, and 
displaying the property for the line. 

Preferably the property relates to nuclear magnetic 
resonances. More preferably, mapping of the property for 
the slice results in a nuclear magnetic resonance image. 

Preferably, the line is a plurality of lines. More 
preferably in the fourth aspect a computer routine maps 
nuclear magnetic resonance values for each of the lines 
by displaying an array of values. 

In one embodiment the object is a biological object such 
as a brain and each of the lines is a line through the 
cortex of the brain. 

Preferably the method further comprises the step of 
analysing the array of values using a fifth computer 
routine. More preferably the fifth computer routine 
analyses numerical properties. 

Preferably analysing the array of values using the fifth 
computer routine comprises differentiation. More 
preferably analysing the array of values comprises 
determining the positions of the maxima, minima and/or 
zero-points of the array of values and the derivatives of 
the array of values. 

In one embodiment, the positions of the inner and outer 
boundaries of the grey matter of the brain may be 
approximated by determining the positions of the maxima 
of the first derivative of the array of values. Two 
closely followed zero-points of the first derivative may 
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indicate the presence of lamination. By measuring the 
relative positions of the zero points the thickness of 
the lamination can be determined. 

According to a sixth aspect of the invention there is 
provided a map produced by a method of mapping a property 
of an object, said method comprising the steps of: 

mapping the property for a slice within the object, 
providing a line within the slice, and 
displaying the property for the line. 



BRIEF DESCRIPTION OF THE DRAWINGS 

Preferred embodiments of the invention will now be 
described, by way of example only, with reference to the 
accompanying drawings in which: 

Figure 1 (a) shows a nuclear magnetic resonance 
image of a living human brain. The insert shows a typical 
line scan. Figure 1 (b) shows lamination within the image 
shown in Figure 1 determined by a computer routine. 

Figure 2 (i) shows a histological image of a section 
of a human brain (intensities changed) , and the insert 
shows the original image. Figure 2 (ii) shows a nuclear 
magnetic resonance image corresponding to the 
histological image of Figure 2 (i) . Figures 2 (iiia) and 
2(iiib) show histological profiles as indicated in Figure 
2 (i) . Figure 2 (iv) shows intensity line scans 
corresponding to regions indicated in Figure 2(i) and 
Figure 2 (ii) . 

Figure 3 shows a mesh of two-dimensional forms that 
approximates the shape of a human brain. 
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Figure 4 shows the analysed nuclear magnetic 
resonance values in a three-dimensional visualisation. 
Darker regions indicate lamination. 

Figure 5 shows a MR I slice of cartilage from a human 
knee and the insert shows line scans along which analyses 
were conducted. 

DETAILED DESCRIPTION OP THE PREFERRED EMBODIMENTS 

Example 1: Two-dimensional Semi -automated Cortical 
Lamination Analysis 

(Computer codes referred to throughout this example are 
listed in the Appendix.) 

High- resolution Tl nuclear magnetic resonance images 
(NMRI) of the brain with high signal, low signal-to- 
noise, and high grey/white and grey/grey matter contrast 
suitable for lamination analysis were obtained. In vivo 
studies were performed using a magnetic field of 1.5 T 
and post-mortem studies were performed at 4.7 T. 

Figure 1 (a) shows a nuclear magnetic resonance image of 
living human brain. Eight Tl weighted spin echo images 
were acquired on a GE Signa NMRI scanner in a single 
session using a 7.5" surface coil positioned at the 
occipital poles. These were aligned, re-sampled and 
averaged. This image shows evidence of the striate cortex 
within the calcarine fissure (arrow) . [The parameters 
were: number of slices = 148; slice thickness = 0.700mm; 
field of view (FOV) = 16 cm 2 ; matrix size = 256 x 256; 
In-plane resolution = 0.625mm x 0.62 5mm; Echo Time (TE) = 
2.7s; Repetition Time (TR) = 12.4s; Number of Excitations 
(NEX) = 1; Flip Angle = 25°; Inversion Time (Tl) = 350ms] . 
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The insert of Figure 2 (i) snows a histological image of 
a section of a human brain and the expanded view shows 
the same image with changed intensities. Figure 2 (ii) 
shows a nuclear magnetic resonance image corresponding to 
the histological image of Figure 2 (i) . [The parameters 
were: number of slices: 30; slice thickness: 1mm (no 
gap); FOV: 8x4 cm2 ; matrix size: 512 x 256; in-plane 
resolution: 156|J,m x 156|xm; TE eff = 82.4 ms; TR = 6 S; NEX 
= 4; Acquisition Time (TA) = 1 hr 42 min 24 s] . Figures 
2 (iiia) and 2(iiib) show histological profiles related 
to areas indicated in Figure 2 (i) . 

Once the images were acquired, they were analysed to 
determine cortical areas of lamination. This involves 
examining the intensity changes across the cortical grey 
matter in the images. It is known that the relative 
concentrations of myelin within cortical laminae can be 
directly measured in the striate cortex using high- 
resolution NMRI . Depending on the scanning protocol used, 
areas of myelination display either a dip or a peak in 
the intensity profile across that area of cortex. 

To the analysis a series of intensity line profiles 
around the entire cortex are taken perpendicular to the 
cortical surface and through the grey matter. Figure 2 
(iv) shows intensity line scans from corresponding 
regions of striate cortex in the histological slide 2 (i) 
and the NMRI slice 2 (ii) and the Insert of Figure 1(a) 
shows a typical line profile showing a characteristic 
intensity drop. ( U A" corresponds to a point just outside 
the cortical surface whilst "B" corresponds to a point 
just deep of the grey/white matter junction) . 
Figure 1 (b) shows lamination determined using the 
computer routine "corla" . This plot shows regions in 
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which there is an intensity dip found within the line 
profile indicating a particular pattern of cortical 
lamination. Results were thresholded and grouped based on 
a measure of the noise from the average image found by 
calculating the standard deviation of the signal 
intensity from an apparently uniform area (the 
background) . Large pale grey circles are those in which 
the intensity drop is between 1 and 2 standard 
deviations, large dark grey circles between 2 and 3, and 
black circles greater than 3 . 

The program "Corla" implements a two-dimensional, slice- 
by- si ice method and requires as input a series of 
contiguous evenly spaced line profiles taken through the 
cortical grey matter, perpendicular to its surface. The 
main steps taken in the analysis consist of determining: 

• cortical boundaries: the air/cortical boundary and the 
grey/white matter boundary, 

• the presence or absence of cortical lamination 
evidenced by an internal intensity drop (ie a local 
intensity peak followed by a trough) - if present, the 
number, position, thickness and intensity drop for the 
lamina are calculated. 

The results are represented as an array of coordinates 
along the cortical boundary each of which has assigned to 
it a series of values which characterise the line profile 
at that point through the grey matter normal to the 
surface. These results are then graphically represented 
for ease of interpretation. 

1 . Line Profiles 

The lines were selected through the cortex and spaced as 
evenly as possibly approximately 3-5 mm apart. Each line 
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was extended passed both cortical boundaries to ensure 
the entire grey matter was sampled. Using the line 
profile function in MEDx3.2, the line-profiles were 
recorded to file. The program to analyse these 
lineprofiles was written in MAT LAB 6.0. The line profiles 
were read in as an array consisting of information about 
each line profile , namely the starting and finishing 
voxels, followed by a series of intensity readings. 

2 . Boundary Determination 

This array is fed into a sub-program, "inflex", which 
determines the inner and outer grey matter boundaries. 
The boundaries are found by assuming they occur where 
there is a maximum change in intensity values between two 
adjacent points ie inflection points in the line profile 
or peaks in the first differential of this plot. There 
will an inflection point at the start and end of each 
line profile marking each boundary. An approximate first 
differential is taken of the line profile by calculating 
differences in intensity values between immediately 
adjacent points. The boundaries are related to the 
maximum peaks in this differential at the start and end 
of the plot. This function finds local maxima by using a 
moving window to smooth smaller fluctuations finds the 
midpoint of any plateaus. The size of the window 
determines the number of distinct local maxima it finds. 
The maximum peak in the first half of the plot is then 
set as one boundary and the maximum peak in the second 
half of the plot as the other. These positions, and their 
corresponding intensity values for all the line profiles 
are saved as an array and are the output of the sub- 
program "inflex" . 
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3 . Lamina Determination 

The original array is also fed into a sub-program, 
"lamf inder" , which determines if any lamination is 
present, the number of lamina, and various 
characteristics of the lamination. Lamination is detected 
by a dip in the intensity profile ie a local peak 
followed by trough. "Lamf inder" finds the peaks and 
troughs by finding where the approximated 1st 
differential crosses zero and where there is a change in 
sign on either side of this point in the first 
differential. Again, some smoothing is done to avoid 
minor non- significant stationary points and to find 
midpoints of plateaus. Once it has determined the 
position of the peaks and troughs and their corresponding 
intensity values, it then calculates: 

• position of the peak and trough relative to the 
cortical boundaries as determined from "inflex". 

• A measure of the thickness of the lamination (ie 
difference in position between peak and trough) 
relative to the cortical thickness. 

• The intensity drop across the lamination 
(intensity difference between peak and trough) 
relative to the intensity drop across the entire 
grey matter (alternatively, relative to the 
intensity value at the outer cortical boundary) 

• Using the output from "inflex", the sub-program 
"coord" calculates the actual (x,y) coordinates of 
the cortical boundaries. These are used for 
plotting the cortical edge to give graphical 
representation of the results. 
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4 . Final Interpretation ana Representation 

Because even minor dips in intensity are detected by this 
program, it is necessary to threshold the results before 
plotting them. The intensity drops across the lamination 
have been used for this purpose. Discarding any 
"laminations" where the intensity drop is less than 5% of 
the total grey matter intensity drop clears a lot of the 
"noise" out of the results, and specific regions of 
lamination begin to appear. Higher thresholds yield more 
localised regions. A number of different aspects of the 
lamination can be plotted. Over a simple 2D plot of each 
slice showing one of the cortical edges the position of 
points of lamination can be plotted. Each point is colour 
coded to represent the degree of intensity drop found for 
that lineprofile. Multiple single slice 2D plots can then 
be stacked to give a 3D plot showing regions of 
lamination on the cortex. Regions with similar intensity 
drops can then be identified, as can a gradual change in 
the intensity drop which may mark a region with poorly 
defined borders. 

Alternatively, using one of the cortical edges as a 
template, a 3D plot can be generated by projecting up 
from this the position of the start and finish of 
lamination found for that lineprofile at that point. This 
is useful in determining how the position and thickness 
of the lamination various around a single slice. 

Automatic line profile generation and analysis 

The following will describe an alternative procedure for 
the generation and analysis of line profiles. This 
process uses a software rountine which can broadly be 
divided into two parts: (A) Automatic line profile 
generation. The operation of the software is now outlined 
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in the context of an analysis of a brain. Line scans that 
were generated and along which analyses were conducted 
using this software routine are shown in Figure 5. 

A. Automatic line profile generation (functions: 
Getbounds; skeleton; skelcontours ; closestpt) 

- calculate a skeletal representation of the grey matter 

of a brain (GM) from a segmented image (use distance 
transform method) 

- the skeleton allows a set of non-intersecting profiles 

to be drawn across the GM 

- these profiles are unevenly spaced and often have many 

share a single end point 

- prune these profiles to keep only the shortest whenever 

they share an end point 

- interpolate between the pruned profiles to get an even 

set of non- intersecting profiles 

Function ; GetBounds 

- Start with the segmented image, with GM=1 # WM=2, 
CSF/other=0 

- calls on: 

Function : skeleton 

- Form a new image such that 

(A) all GM points that border the with matter of a 
brain (WM) and let them=2 

(B) also find all GM points that border CSF/other 

and let them=3 

(C) let all other GM points=l 

(D) and all remaining (non-GM) points=0 



Func t i on : skelcontours 
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- For all points that=l in this new image find 

(using function: closestpt) 

(A) the minimum distance for a point of value 2 

(B) the minimum distance for a point of value 3 

- Whenever the above two distances differ by a small 

amount (say sqrt(2) in 2D) then that point is 
marked as a skeleton point and the points of 
values 2 and 3 that are at the minimal 
distances are the end points of the profile. 
Prune this set by only keeping the shortest 
profile 

out of each set of profiles that share an end 
point . 



B. Line profile analysis (Functions : TwoD; iprofile; IPA; 
coord; stationary/ concavity; peaks) 

Function; TwoD 

- read in the coordinates of the line profiles generated 
by Getbounds 

- user generated input read in 

- calls on functions: 

Function: iprofile 

generates interpolated line profiles between the 
specified coordinates 

Function: IPA 

returns the coordinates of the stationary points of 
the line profile / widths, relative depth and 
intensity drops related to the stationary points 
(using function: stationary; peaks; coord) 
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returns the coordinates of the concavity changes of 
the line profile (using function: concavity; peaks, 
coord) 

plots the line profiles if required 

- Output is user specified results overlaid on original 
image using data returned: choices of results to view 
include position, depth, number, width and intensity 
drop of laminations. 



Example 2: Three-dimensional Volumetric Automated 
Cortical Lamination Analysis 

[Computer codes referred to throughout this example are 
listed in the Appendix. An off the shelf function 
"spm_hread" was used to read in the header information 
for the image (SPM99, freely available) . All computer 
codes are in MATLAB 6.0] . 

The 3D analysis builds on the 2D analysis already 
described. In this example there are five steps in the 3D 
approach . 

1. 3-D Representation of the Brain 

As line profiles need to be dropped perpendicularly 
through the cortex, the surface of the brain needs to 
be accurately characterised in order to get the correct 
topography. 

The surface was represented as a mesh, using the 
software package EMSE Suite (Source Signed Imaging, 
Inc.), although the analysis is not dependent on what 
program was used to generate the mesh, and can be 
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adapted to use others. EMSE approximates the surface as 
a mesh made up of triangles. 

2. Information (meshes and images) are read in and normals 
to surfaces calculated 

This module needs to read EMSE (or in the future 
other formats) meshes. It outputs lines along which the 
profiles need to be generated. Initially it has been 
agreed that the normals to the faces of the mesh would 
be calculated. This is due to coding and computational 
simplicity. The degree of sampling therefore depends on 
the size of the mesh. As the number of faces is 
increased, their size decreases and hence the density 
of the sampling will also increase. Additional normals 
at vertices, and further interpolation along the edges 
and/or faces can then be added incrementally so that 
ultimately, every point on the surface of the cortex 
can be sample. "Quickrun" is the initiating program, 
prompting the user to input inf ormation and choose 
analysis preferences. It calls the function "read_emse" 
which reads in the mesh. The functions "tricenter" , 
"trinormal" and "normalize" then calculate the normals 
to each triangular face of the mesh through its centre. 
The function "snc" is then called which generates the 
line profiles and performs the analysis. 

3 . Use normals and original Analyze image to calculate 
intensity line profiles 

This is the function of "snc" . It reads in the image 
data using "anzload" and, using the normals calculated 
previously, generates the line profiles. The function 
"Ipa" is then called which performs the actual line 
profile analysis. 
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4. Perform cortical lamination analysis on line profiles 

u Lpa" determines changes in concavity stationary 
points and boundary points by calling various 
functions. The function "concavity" finds the changes 
in concavity along the line profile. These correspond 
to where the second differential crosses zero. 
Therefore, it passes the first differential of the line 
profile to the function "peaks" which then finds the 
differential of this (ie it uses the second 
differential in this case) . The zero points are found 
by determining where there is a change in sign along 
this curve. The function ""boundary_j?ts" then uses the 
information from "concavity" to determine the cortical 
boundaries using a vector method which incorporates a 
number of competing criteria for assigning the cortical 
boundaries. The function "stationary" determines 
stationary points in the line profile. These correspond 
to where the first differential crosses zero. 
Therefore, "stationary" passes the line profile to 
"peaks" which then finds the zero points in the first 
differential in the same way it did for "concavity" . 
Using all this information, "lpa" then determines a 
number of features of each line profile such as the 
number of laminations, their widths and the intensity 
drops corresponding to each. 

All data is thresholded using an intensity change 
significance threshold and a resolution significance 
threshold. The former is determined from the signal-to- 
noise ratio of the original image whilst the latter is 
determined from the resolution of the original image. 



5. Visualisation of results 

The information from "lpa" is used by "snc" to 
colour the mesh using the function "colour_picker" 
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according to what results the user wishes to visualize. 
A colour coded wireframe in OOGL (freeware: suitable 
for geomview) is generated. OOGL can be easily- 
converted to VRML using available tools. This produces 
a 3D model which can be rotated and rendered in real 
time using geomview (freeware) . The colour represents a 
particular analysis result and the intensity of that 
colour quantitatively reflects the magnitude of that 
result. Fanother option may involve overlaying the MR 
and the wireframe and also combining multiple coloured 
meshes representing various combinations of results. 
Similarly coloured regions would represent cortical 
regions sharing a common lamination feature eg, 
intensity drops, positions of laminae, thicknesses, 
number of laminae etc. Figure 3 shows a mesh of two- 
dimensional forms that approximates the human brain. 
Darker regions correspond to darker shades of colour 
and indicate, in this example, lamination. Figure 4 
shows the corresponding three-dimensional visualisation 
without the mesh. 

It will be appreciated that each line profile contains a 
large amount of information about the cortex through 
which it passes in the MR I image. In this example it has 
been chosen to examine one aspect of that by using first 
and second differentials to determine stationary points 
and changes in concavity. As the quality and or type of 
images improves, or even regardless of this, there are 
many other ways of examining and analysing the 
information that can be incorporated into this analysis 
which can aid in defining distinct regions. 

It will also be appreciated that this method is not 
limited to biological objects. The method of mapping the 
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property for the object is not limited to the detection 
of a mapping the property for a plurality of slices. 
Alternative methods can be used which may not comprise 
the detection of slices. It will also be appreciated that 
the method is not limited to nuclear magnetic resonance 
imaging. Each of the lines may also be a region. 
Alternative computer routines may be used to analyze the 
array of values which may or may not comprise numerical 
procedures such as differentiation. 

It will also be appreciated that a procedure similar to 
the automatic line profile generation and analysis 
described in the context of 2-D analysis can be extended 
to be applicable 3-D analysis. 

Although the invention has been described with reference 
to particular examples, it will be appreciated by those 
skilled in the art that the invention may be embodied in 
many other forms . 

It is to be understood that, if any prior art publication 
is referred to herein, such reference does not constitute 
an admission that the publication forms a part of the 
common general knowledge in the art, in Australia or any 
other country. 
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Appendix 

Codes of computer routines referred to throughout the 
examples 

Computer code used for 2-D analysis 

function [I^m3 5 LamData,Idrop]==Corla(tl,t2,fig,Z) 

[lpStart,lpFinish,x]=reader3; 

[r,n]=size(x); 

y=([i:rD; 

LamData=[]; 
fori=l:n 

z=[y,x(:,i)]; 

t=diff(z(:,2)); 

[E(i).ce,E(i)gwb]=inflex(z); 
[Lam(i).peak,Lam(i).tro^ 
if isempty(LamD) 

LamData(:,:,i)=zeros(30,3); 

Idrop(:,:,i)=zeros(30,3); 
else 

iflength(LamD(:,l)>— 30 
Id(30,3)=0; 

LamD(30,3)=0; 
end 

LamD ata( : , : , i)=LamD ; 
Idrop(:,:,i)=Id; 
end 

Profileplotter(z,E(i),Lam(i),i); 
end 

[Start ,Finish]^oord(LamData,lpStart,lpFinish,E,tl,t2,fig,Idrop,Z); 



function [Start,Finish]=coord(LaniData4pStart,lpFinish,E,tl,t2,fig,Idrop,Z) 

Start=[]; 

Finish=[]; 

figure(fig); 

hold on 

for i=l :length(lpStart) 

[S,F]=pointplot(lpStart(i,:),lpFinish(i,:),E(i)); 
Start=[Start;S]; 
Finish=[Finish;F]; 
plotCStartCU^StartCUXW); 
hold on 
end 

h^findobjCgc^TypeMineO; 
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set(hl(:), , MarkerSize , ,3.5); 
setOilCO/LineWidth'^.S); 
for i=l :length(lpStart) 
if (~isempty(find(Idrop(: ,3,i)>=t 1 ))) 

plotCStartCUXStartCUVro'); 

hold on 

h^frndobjCgcf/Typeyiine'); 
set(hl(l), , LineWidth , ,7.5); 
end 
end 

for i=l :length(lpStart) 
if (~isempty(find(Idrop(:,3,i)>=t 1 ))) 
plotCStartCUXStartO^Vgo'); 
hold on 

h^findobjXgc^Type'/line 1 ); 

set(hl(lVMarkerSize\5.5); 
set(hl(lVLineWidth\5.5); 
end 
end 

for i=l :length(lpStart) 
ifHsemptyCfindCIdropf^sa^^^tl)))) 
plot(Start(i, l),Start(i,2), , yo'); 
hold on 

hl=findobj(gcf,Typeyiine'); 

set(h 1 ( 1 ),'MarkerSize\5 . 5) ; 
set(hl(lVLineWidth',5.5); 
end 
end 

for i=l :length(lpStart) 
ifHsempty(find(Idrop(:,3,i)>=(3*tl)))) 
plot(Start(i,l) > Start(i,2), , ro'); 
hold on 

hl^findobjCgc^Type'/line 1 ); 

set(hl(l)/MarkerSize\5.5); 
set(hl(iyLineWidth\5.5); 
end 
end 

figure(fig+l); 
hold on 

for i=l:length(lpStart) 
t=max((Idrop(:,3,i))); 
v=find(Idrop(:,3,i)=t); 
ift>3 

plot3(Start(i,l),Start(i,2),LamData(v(l),2,i), f b*'); 
hold on 

ploG(Start(i,l),Start(i,2),LamData(v(l),l,i), , r* , ); 
hold on 
end 
end 
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%g=i; 

% while (g<=length(LamData(:,l,i)))&(LamData(g,l,i)>0) 

% if (LamData(g,2,i)-LamData(g, 1 ,i))>5 

% plot3(StaIt(i,l),Start(i,2),I^mData(g,2,i), l b* , ); 

% hold on 
%plot3(Start(i, 1 ),Start(i,2),LamData(g, 1 ,0,'**'); 

% hold on 

% end 

% g=g+l; 

% end 
%else 

% plot3(Start(U),Start(i,2),Z, W); 
% hold on 



function [S,F]=pointplot(lpStart,lpFinish,E) 

S(l^)-lpStart(l,l)+(Exe(l)7200)*(lpFinish(l,l)-lpStart(l,l)); 
S(l,2)=lpStart(l,2)+(E.ce(l)/200)*(lpFinish(l,2)-lpStart(l,2)); 
F(l,l)=lpStart(l,l)+(E.gwb(l)/200)*(lpFinish(ia)-lpStart(l,l^ 
F(l,2)=lpStart(l,2)+(E.gwb(l)/200)*(lpFinish(l,2)-lpStart(l,2)); 



function [PPos,Pheight]=avg(lp,x,PPos) 
tsum=0; 

for i=l:length(x) 
ifPPos<8 

tsum=tsum+x(i); 
else 

tsum=tsum+(x(i)-9) ; 
end 
end 

PPos=round(((PPos*length(x))+tsum)/length(x)); 
Pheight=lp(PPos); 



function [ce,gwb]=inflex(lp) 

x=diff(lp(:,2)); 

n=length(x); 

P=[l,x(l)]; 

count=0; 

i=l; 

d=l; 

inf=[]; 



\ 
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HI; 

while i<n-10 
bit=x(i:i+9); 
localmax=max(bit); 

y=find(bit>(localmax- 1 )); 

temp=size(y); 

if count<=l 

d=y(temp(l,2)); 

else 

d=y(D; 

end 

ifx(d+i)>P(l,2) 

P=[d+i,x(d+i)]; 

else 

if count==0 

P=refiner(x,P,n); 

inf=[inf;P]; 
count=count+l; 
elseif (d+i)-inf(count, 1 )> 1 0 
P=refiner(x,P,n) ; 
inf=[inf;P]; 
count=count+l; 
end 

P=[d+i,x(d+i)]; 

end 
i=i+l; 

end 

a=find(inf(:,l)<100); 
[ceI,ceP]=max(inf(l :length(a),2)); 
ce=[infl[ceP, 1 ),cel] ; 

[gwbI,gwbP]=max(inf(length(a)+l:coiint,2)); 
gwb=[inf((gwbP+length(a)) 5 1 ),gwbl] ; 
I=[lp(inf(:,l)),inf][:,2)]; 



function [Peak,Trough,LamData,Idrop]=lamfinder(lp,start,finish) 
if finish<length(lp); 

n=finish; 
else 

n=length(lp); 
end 

count=0; 
temp=0; 
i=20; 

f=finish; 

x=diff(lp); 

peakp=[]; 
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peakl=[]; 
troughp=[]; 
troughI=[]; 
while i<(n(l)-31) 
if(x(i))<=0 
temp=i; 

while (x(temp)=0)&(temp<(n(l)-31)) 

terap=temp+l; 
end 

ifx(temp)<0 

count=count+l; 
peakp=[peakp ;(temp+i)/2] ; 

peakl=[peakl;lp(temp)] ; 
while (x(temp)<0)&(temp<(n( 1 )-3 1 )) 

temp=temp+l; 

end 

i=temp; 

if x(temp)=0 

while (x(temp)=0)&(temp<(n( 1 )-3 1 )) 

temp=temp+l; 

end 

end 

troughp=[troughp ;(i+temp)/2] ; 
troughI=[troughI;lp(i)] ; 
else 

i=temp+l; 
end 
else 

W+l; 
end 
end 

Peak=[peakp,peakl] ; 

Trough=[troughp,troughI] ; 

Pdata=[]; 

Tdata=[]; 

Idropl=[]; 

Idrop2=[]; 

if ~isempty(Peak) 

Pdata^(Peak(:, l)-start)* 1 00)/(finish-start); 

Tdata={(Trough(: , 1 )-start)* 1 00)/(finish-start); 

LamData=[Pdata,Tdata] ; 

Idrop3=<Peak(:,2)-Trough(:,2)); 

IdroplK(Peak(:,2)-Trough(:,2))*100)/(lp(finish)-lp(start)); 

Idrop2=((Peak(:,2)-Trough(: > 2))*100)/lp(start); 

Idrop=[Idrop 1 ,Idrop2,Idrop3] ; 
else 

LamData=[]; 

Idrop=[]; 
end 
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function ProfileplotterfoEjLanv); 
figure(i); 

grid on 
hold on 

plot((E.ce(l)),(E.ce(2));ro',(E.gwb(l)),(E.gwb(2));ro , ) 
hold on 

if -isempty(Lam.peak) 

plot((Lam.peak(l , 1 )),(Lam.peak( 1 ,2)),'ro'); 

plot((Lam.trough( 1 , 1 )),(Lam.trough( 1 ,2))/ro'); 
end 



function [Start,Finish,mat]=reader 

Start=[]; 

Finish=[]; 

mat=[]; 

h=dir( , 106-7* t ); 
n=length(h); 
for i=l:n 
%s= s h(i).name 

fid=fopen(h(i).name,'rt'); 
forj=l:2 

fgetl(fid); 
end 

crap=fscanf(fid, , %s\t%s\t%s\t, , ); 

SKfscanf(fid/%f,'))'; 

fgetl(fid); 

crap=fscanf(fid, , %s\t%s\t%s\t, 1 ); 
F=(fscanf(fid > , %f;)) , ; 
Start=[Start;S]; 
Finish=[Finish;F] ; 
for j=l:6; 

fgetl(fid); 
end 

temp=[]; 
forj=l:200 

data=fscanf(fid,'%i\t%i,'); 

fgetl(fid); 

temp=[temp;data(2, 1 )] ; 
end 

%temp2=[]; 
%forj=l:200 

% temp2=[temp2;temp(201-j)]; 
%end 
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mat=[raat,temp] ; 
fclose(fid); 
end 



function Peak=refiner(lpJ\n) 

PPos=P(l); 

Pheight=lp(P(l)); 

if (PPos>8)&(PPos+8<n) 

[Pheight,temp]=max(lp(PPos-8:PPos+8)); 

PPos=temp-9+PPos; 
elseifPPos<=8 

[Pheight,PPos]=max(lp(l :PPos+8)); 
else 

[Pheight,temp]=max(lp(PPos-8:n)); 
PPos=temp-9+PPos; 
end 

if (PPos>8)&(PPos+8<n) 

x=find(lp((PPos-8:PPos+8))>Pheight-3); 

[PPos,Pheight]=avg(lp,x,PPos); 
elseifPPos<=8 

x=find(lp((l :PPos+8))>Pheight-3); 

[PPos,Pheight]=avg(lp,x,PPos); 
else 

x=find(lp((PPos-8:n))>Pheight-3); 
[PPos,Pheight]==avg(lp,x,PPos); 
end; 

Peak=[PPos,Pheight]; 



Alternative computercode used for 2 -D analysis 

function [imO, coords] = Getbounds 
nl = sprintf ( ' \n • ) ; 

disp(* Please enter the name of the file to be analysed'); 
Image= deblank ( input ( ' Image file: ','s 1 )); 
disp (nl) ; 

disp(' Please enter the name of the segmented grey matter file 1 ); 
Grey=deblank (input ( 'Grey matter image f ile : ' , ' s ' ) ) ; 
disp (nl) ; 

disp ('Enter 1 if the file is avw, and 0 if it is ppm'); 
ftype = input ('file type: *); 
disp (nl) ; 

disp ('Do you want to subsample the image?'); 

disp ('(type 1 if you do not, otherwise the factor to downsample by 
(integer greater than 1)) '),- 
sf = input (' subsample factor: '); 
disp (nl) ; 
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disp(* Please enter the name of the output file for the 
coordinates? • ) ,* 

base = deblank (input (• file name: ' ,'s')); 
disp (nl) ; 

if ftype==0, 

imO = readppm ( Image ) ; 

gm = readppm (Grey) ; 
elseif f type==l , 

imO =read_avw( Image ) ; % read original image 

gm=read_avw (Grey) ; %read grey matter image 

end 

imO = imO ( : , : , 1) ; 
gm = gm ( : , : , 1 ) ; 
figure (1) ; 
image sc (imO) ; 
colormap(gray) ,- 
hold on; 

% now make the segmented image (WM=2, GM=l, other=0) 
segim = (im0>l) *2 - (gm>l) ; 

% subsample the image (factor of "sf") 
if sf > 1 

[N,M] =size (segim) ; 

segim = segim (l:sf:N,l:sf:M); 

end 

figure (2) 
image sc (gm) ; 
colormap (gray) ; 
pause 

figure (3) 
imagesc (segim) 
colormap (gray) ; 
pause 

% fix up any problems (I had to fill a couple of holes in the GM) 
% now generate the coordinates of the unique profiles 
coords = getcontours (segim) ; 

% convert coordinates back to original if necessary (upsample) 
if sf > 1 

coords = (coords - 1) *sf + 1; 

end 

% write the coordinates to file 
output^ file = strcat(base, '.dat'); 
fout = f open (output_f ile, 1 w 1 ) ; 
for i=l : length (coords) 

fprintf (f out, 1 %f %f %f %f\n', coords (i,l), coords (i, 2), 
coords (i, 3), coords (i , 4 )) ; 

end 

f close (fout) ; 



function [dims , scales , bpp, endian, datatype] = read_avw_hdr (f name) 

% [dims , scales , bpp , endian] = READ_AVW_HDR ( f name ) 

% 
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% Extracts the 4 dimensions (dims) , 

% 4 scales (scales) and bytes per pixel (bpp) for voxels 
% contained in the Analyse header file (fname) 

% Also returns endian = ' 1' for little-endian or , b» for big-endian 
% 

% See also: READ_AVW, READ_AVW_IMG, SAVE_AVW, S AVE_AVW_HDR , 
SAVE_AVW_IMG 

% remove extension if it exists 
if ( (length (findstr (fname, * .hdr' )) >0) | ... 
(length(f indstr ( fname, ' .img' ) >0) ) ) , 
fname=f name (1 : (length (fname) -4) ) ; 
end 

fnhdr=s treat (fname, ' .hdr 1 ) ; 

% open file in big-endian 
endian= *b » ; 

f id=fopen(fnhdr, 'r • , 'b' ) ; 
testval = f read (fid, 1, ' int32 • ) ; 

% check if this gives the correct header size - if not use little- 
endian 

if (testval-=348) , 
f close (fid) ; 

f id=fopen(fnhdr, 'r • , ' 1 1 ) ; 
endian= 1 1 ' ; 

testval = fread(fid, 1, 'int32 ' ) ; 
if (testval~=348) , 

disp('Can not. read this file format'); 

return; 
end 
end 

% ditch the remaining initial header stuff 

dummy =f read(f id, 36, 'char') ; 

% ditch dim[0] = No. dimensions 

dummy=f read (fid, 1, 1 intl6 ' ) ; 

dims=f read (fid, 4 , * intl6 ' ) ; 

dummy=fread(fid,3, 'intl6') ; 

dummy=f read (fid, 14 , ' char ' ) ; 

datatype=f read (fid, 1, f intl6») ; 

bpp=f read (fid, 1, ■ intl6 « ) ; 

dummy^f read (f id, 2 , 'char'); 

dummy=f read (fid, 1, 'float 1 ); 

scales=f read (fid, 4 , • float ! ) ; 
f close (fid) ; 
return; 



function img = read_avw_img (filename) ; 

% [img] = READ_AVW__IMG (filename) 

% 

% Read in an analyse file into either a 3D or 4D 

% array (depending on the header information) 

% Ouput coordinates are in MEDx convention 

% except that all dimensions start at 1 rather than 0 

% Note: automatically detects char, short, long or double formats 

% 

% See also: READ__AVW , READ_AVW_HDR, SAVE_AVW, S AVE_AVW__HDR , 
SAVE AVW IMG 



% remove extension if it exists 

if ( (length (findstr (filename, » .hdr' )) >0) | 
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(length (findstr (filename, ' . img • ) >0) ) ) , 
filename** filename (1 : (length (filename) -4) ) ; 
end 

fnimg=s treat (filename, 1 . img' ) ; 
fnhdr=strcat (filename, 1 -hdr') ; 

[dims , scales, bpp,endian, datatype] = read_avw_hdr (f nhdr) ; 
fp=f open (f nimg, 'r' , endian) ; 
if (datatype==4 ) , 

dat=f read (f p, ' short • ) ; 
elseif (datatype==2) , 

dat=f read (fp, ' char 1 ) ; 
elseif (datatype==8) , 

dat=fread(fp, ' int • ) ; 
elseif (datatype ==64 ) , 

dat =f read (fp, 'double 1 ) ; 
elseif (datatype==16) , 

dat=fread(fp, 'float32 ' ) ; 

end 

f close (fp) ; 

nvox = prod (dims) ; 

if ( length (dat) <nvox) , 

error (' Cannot open image as .img file does not contain as many 
voxels as the .hdr specifies'); 
elseif (length (dat) >nvox) , 

disp ( 'WARNING: : truncating .img data as it contains more voxels than 
specified in the .hdr'); 

dat - dat (1: nvox); 
end 

if (dims (4) >1) , 

img = reshape (dat, dims (1) , dims (2) , dims (3) , dims (4) ) ; 
else 

img m reshape (dat , dims (1) , dims (2) , dims (3) ) ; 
end 

clear dat; 

%% DEFUNCT FLIPPING 

%% flip y dimension to be consistent with MEDx 
%img=f lipdim(img, 2) ; 



function [coords] =getcontours (segim) 

% coords = getcontours (segim) 

% 

% Input (segim) is a segmented image (2D) which has 

% values of 1 for GM, 2 for WM and 0 otherwise 

% Output (coords) is a matrix where each row contains the 

% end coordinates of a cortical profile (4 values) 

% The intensity profile can be extracted using 

% the iprofileO function with a given row from coords 

% 

% See also IPROFILE, SKELETON, SKELCONTOURS 

im=conv2 (segim, ones (3 , 3) /9, 1 same 1 ) ; 
im2= (im<0 . 9) . * (segim==l) ; 
im3= (im>l.l) . * (segim==l) ; 
im= (segim==l) +2*im2+im3 ; 
sk=skeleton (im) ; 
coords=skelcontours (im, sk) ; 
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function [sk] =skeleton (im) 

% sk=skeleton (im) 

% 

% Calculates the skeleton of an image im (medial axis transform) 
% The image im is a binary image of the area with voxels at 
% one boundary set to 2 and at the opposite boundary set to 3 
% e.g. GM/WM interface = 2, GM/CSF interface = 3 , GM = 1 

[N M] =size (im) ; 
skd=zeros (N, M, 2) ; 

skd(:,:,l) = im*0 + max (M A 2 , N A 2 ) +1 ; 
skd(:,:,2) = im*0 + max(M A 2 / N A 2) +1; 

% Calculate the minimum distance transform 
xim=kron( (1:N) • , ones (1,M) ) ; 
yim=kron(l :M, ones (N, 1) ) ; 
val=[2 3] ; 
for r=l:2, 
for x=l:N, 
for y=l:M, 

if (im(x,y) ==val (r) ) , 

dist = (xim-xim(x,y) ) . A 2 + (yim-yim (x,y) ) . A 2 ; 
skd(: / :,r) = min(skd( : , : , r) ,dist) ; 
end 
end 

disp(x) 
end 
end 

sk = (abs (sqrt (skd( :,:,!)) -sqrt (skd( : , : , 2) ) ) <sqrt (2) ) . * (im==l) ; 



function [coords] =skelcontours (im, sk, cand) 

% [coords] =skelcontours (im, sk) 

% 

% Calculates the *unique* line profile contours from a skeleton (sk) 
% of an image im (medial axis transform) - as obtained from 
% skeleton () 

% Returns the coordinates as an N by 4 matrix - each row having 
% the form [xl yl x2 y2] 

% The image im is a binary image of the area with voxels at 
% one boundary set to 2 and at the opposite boundary set to 3 
% e.g. GM/WM interface = 2, GM/CSF interface =3, GM » 1 
% 

% See also SKELETON 

[sx,sy] = find(sk==l); 
N = length (sx) ; 

% get candidate coord pairs 
if (nargin<3) , 

cand= [ ] ; 

im2= (im==2) ; 

im3= (im=-3) ; 

disp(N) ; 

for k=l:N, 

cand=[ cand ; closestpt (im2 , [sx (k) sy(k)]) closestpt (im3 , [sx (k) 
sy(k)]) ]; 

disp(k) ; 
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end 
end 

% eliminate multiple matches 

for r=l : 2 , % r=l => boundary with im==2 , r=2 => bdy with im==3 
coords = t] ; 

for k-1: length (cand( : , 1) ) , 

if (min(cand(k, :))>0), % only do valid end points 

diff=abs(cand(: , (r*2-l) )-cand(k, (r*2-l) ) ) + abs(cand(: ,r*2) - 
cand(k,r*2) ) ; 

idx=f ind(diff ==0) ; % all rows with same end point 
mindiff = N*2; 
minidx = 0 ; 
mincand= [0 0 0 0] ; 

% find minimum contour with this end point 
for m=l : length (idx) , 

d - norm( [cand(idx(m) ,3) -cand(idx(m) , 1) ,cand(idx(m) , 4) - 
cand(idx(m) ,2) ] ) ; 

if (mindif f >d) , 

mindiff = d; 

minidx = idx(m) ; 

mincand = cand (minidx, :) ; 

end 

cand(idx(m) , : ) = [0 0 0 0]; 
end 

if (minidx>0) , 

if (min (mincand) >0) , 

coords= [coords ; mincand] ; 

end 
else 

error ( 1 FAILED ASSERTION: minidx>0 ' ) ; 
end 

disp(k) ; 
end 
end 

cand = coords; 
end 



function [xc] =closestpt (im,x) 

% xc=closestpt (im,x) 

% 

% Calculates the closest voxel coordinate to x where the 

% binary image is 1, (i.e. arg min_{xc} norm(x-xc) st im(xc)=l) 

[N M] =size (im) ; 

% Calculate the minimum distance transform 

maxd= (M+N+l) *2; 

xim=kron( (1:N) 1 , ones(l,M) ) ; 

yim=kron (1 :M, ones (N, 1) ) ; 

dist = im.* ( (xim-x(l) ) . A 2 + (yim-x (2) ) . A 2) + maxd* (1-im) ; 
[cx, cy] =f ind (dist==min3 (dist) ) ; 
xc(l)=cx(l) ; 
xc(2) sC y(l) ; 



%%%%%%%%%%%%%%%%% 
function TwoD(imO) 
%%%%%%%%%%%%%%%%% 
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nl = sprintf ( ' \n' ) ; 

disp ( • Please enter the name of a coords file'); 
fin = deblank ( input ( 'coords file: ' ,'s')) 
disp(nl) ; 

disp (' Standard deviation of background noise?'); 

disp(' (Enter the value 1, to ignore signifigance rel . to background 
noise. ) • ) ; 

std_dev = input ('std. deviation: '); 
disp (nl) ; 

disp ('Enter value for number of points to take per voxel'); 
pspacing = input ( * spacing : • ) ; 
disp (nl) ; 

disp ('Do you want to display the line profiles?'); 

disp(' (type 0 if you do not, any other number if you do.) '); 

pi = input ( 'display profiles: '); 

disp (nl) ; 

disp ('If you want to plot the number of laminae, enter 0'); 
disp ('If you want to plot the position and intensity drop of the 
myelinated laminae, enter 1'); 

disp ('If you want to plot a combination of number of laminae with the 

position, and intensity drops of the myelinated laminae, enter 2 1 ); 

disp ('If you want to plot the relative centre of each myelinated 

laminae, enter 3'); 

result = input ( 'option: 1 ) ; 

disp (nl) 



fid = f open (f in); 
line « fgetl (fid) ; 
C_num=0 ; 

while line ~= -1, 

C = sscanfdine, ' %f %f %f %f » ) » ; 
line = fgetl (fid) ; 
C__num=C_num+l ; 
coords (C_num, : ) = [C] ; 

end 

f close (fid) ; 

[t , u] =size (coords) ; 

figure (10) ; 

image sc (imO) ; 

colormap (gray) ; 

hold on; 



if result==0 
for i= l:t 

ip = iprof ile (coords (i, : ) , imO, pspacing) ; 
[Ledges, Lcentres, widths, num_ lam, Id,W,relW]= 
IPA (coords (i, :), ip* , std_dev, pi, pspacing, (i+1) ) ; 
figure (1) 
if nutn_lam <=1 

plot ( [coords (i, 2) coords (i, 4) 3 , [coords (i, 1) 
coords (i, 3) ] , ' r-o' ) ; 

elseif num_l am <=2 

plot ( [coords ( i , 2 ) coords ( i , 4 ) 3 , [coords ( i , 1 ) 
coords (i, 3) 1 , 'y-o ' ) ; 

elseif num_lam <=4 

plot ( [coords ( i , 2 ) coords ( i , 4 ) ] , [ coords ( i , 1 ) 
coords (i, 3) 3 , 'g-o' ) ; 

elseif num_lam <=6 

plot ( [coords (i , 2) coords (i, 4) ] , [coords (i, 1) 
coords (i, 3) 3 , 'c-o' ) ; 
else 
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plot ( [coords (i, 2) coords (i # 4) ] , [coords (i, 1) 
coords (i, 3) ] , 'b-o ' ) ,- 
end 

hold on; 

end 

elseif result == 1 
for i= l:t 

ip = iprofile (coords (i, :), imO , pspacing) ; 
[Ledges, Lcentres, num_lam, Id, pcount , W, relW] = 
IPA (coords (i, :), ip ' , std_dev, pi, pspacing, (i+3) ) ; 
figure (1) 

plot ( [coords (i, 2) coords (i, 4) ] , [coords (i, 1) coords (i , 3 ) 3 , ■ r- ' ) ; 
hold on; 

for j=l: pcount 

if Id(j,2) < 1 

plot (Lcentres (Id(j , 1) ,2) , Lcentres (Id ( j , 1) ,1) , , r+' 
elseif Id(j,2) < 2 

plot (Lcentres ( Id (j,l) ,2) , Lcentres (Id (j , 1) , 1) , 'y+' 
elseif Id(j,2) < 3 

plot (Lcentres (Id (j,l) ,2) , Lcentres (Id (j , 1) , 1) , 'y*' 
elseif Id(j,2) < 4 

plot (Lcentres (Id (j , 1) , 2) , Lcentres ( Id ( j , 1) , 1) , ' g+ 1 
elseif Id(j,2) < 5 

plot (Lcentres (Id (j,l) ,2) , Lcentres (Id (j , 1) , 1) , 'g*' 
elseif Id(j,2) < 6 

plot (Lcentres (Id( j , 1) , 2) , Lcentres ( Id ( j , 1) , 1) , ' c+ ' 
elseif Id(j,2) < 7 

plot (Lcentres (Id(j , 1) ,2) , Lcentres (Id (j , 1) , 1) , 'c*' 
elseif Xd(j,2) < 8 

plot (Lcentres (Id(j , 1) ,2) , Lcentres ( Id (j , 1) ,1) , »b+' 
elseif Id(j,2) < 10 

plot (Lcentres (Id(j , 1) ,2) , Lcentres ( Id (j , 1) , 1) , 'b*« 

else 

plot (Lcentres (Id(j , 1) ,2) , Lcentres (Id (j , 1) ,1) , , m+' 

hold on; 
end 

end 

end 

elseif result==2 
for i= l:t 

ip = iprofile (coords (i, :), imO, pspacing) ; 
[Ledges, Lcentres, num_lam, Id, pcount , W, relW] = 
IPA(coords (i, : ) , ip 1 , std_dev, pi, pspacing, (i+3) ) ; 
figure (1) 
if num_lam <=1 

plot ( [coords (i, 2) coords (i,4) ] , [coords (i, 1) 
coords (i, 3) ] , 'b- • ) ; 

elseif num_lam <=2 

plot ( [coords ( i , 2 ) coords ( i , 4 ) ] , [ coords ( i , 1 ) 
coords (i, 3) ] , 'g- ' ) ; 

elseif num_lam <=4 

plot ( [coords (i, 2) coords (i, 4) ] , [coords (i, 1) 
coords ( i , 3 ) ] , 1 c - 1 ) ; 

elseif num_lam > 4 

plot ( [coords (i, 2) coords (i # 4) ] , [coords (i, 1) 
coords (i, 3) ] , 'w- . ■ ) ; 
end 
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hold on; 

if num_lam <=4 

for j=l:pcount 

if Id(j,2) < 4 & Id(j,2) > 1 

plot(Lcentres(W(j / l) ,2) , Lcentres (W ( j , 1) , 1) , ' r* ' ) ; 

hold on; 

plot (Lcentres (W(j , 1) , 2) , Lcentres (W ( j , 1) , 1) , ■ ro 1 ) ; 
elseif Id(j,2) < 7 

plot (Lcentres (W(j f l) ( 2) , Lcentres (W(j , 1) , 1) , 'y* 1 > ; 
hold on; 

plot (Lcentres (W(j,l) ,2) , Lcentres (W (j , l) f 1) f , yo , > ; 
elseif Id(j,2) < 15 

plot (Lcentres (W(j # l) ,2) , Lcentres (W (j , 1) ,1) , 'w*') ; 
hold on; 

plot (Lcentres (W( j , 1) , 2) , Lcentres (W (j , 1) , 1) , 'wo' ) ; 
elseif Id(j,2) >=15 

plot (Lcentres (W( j , 1) , 2) , Lcentres (W ( j , 1) , 1) , 1 w* • ) ; 
hold on; 

plot (Lcentres (W(j , 1) ,2) , Lcentres (W (j , 1) ,1) , *ws') ; 

end 
hold on; 
end 
end 

end 

elseif result==3 
for i= l:t 

ip = iprofile (coords (i, :), imO ,pspacing) ; 
[Ledges, Lcentres, num_lam, Id,pcount, W, relW] = 
IPA (coords (i, :), ip* , std_dev, pi, pspacing, (i+3 ) ) ; 
figure (10) 

plot ( [coords (i, 2) coords (i,4) ] , [coords (i , 1) coords (i , 3 ) ] , 1 r- • ) ; 

hold on; 

if num_lam <=4 

for j=l:pcount 

if relW(j,2) <= 0.15 

plot (Lcentres (relw( j , 1) , 2) , Lcentres (reiw ( j , 1) , 1) , 'm* ' ) ; 

plot (Lcentres (relW( j , 1) , 2) , Lcentres (relW( j , 1) , 1) , 1 mo ' ) ; 
elseif relW(j,2) <= 0.325 

plot (Lcentres (relW(j,l) ,2) , Lcentres (relW (j , 1) , 1) , • b* ' ) ; 

plot (Lcentres (relW(j,l) ,2) , Lcentres (reiw (j , 1) ,1) , 'bo') ; 

elseif relW(j,2) <= 0.5 

plot (Lcentres (relW( j , 1) , 2) , Lcentres (relW (j , 1) , 1) , ■ c* * ) ; 

plot (Lcentres (relW(j , l) ,2) , Lcentres (relW (j , 1) , 1) , 'co'); 
elseif relW(j,2) <= 0.675 

plot (Lcentres (relW(j , 1) ,2) , Lcentres (relW (j , 1) ,1) , 'g*') ; 

plot (Lcentres (relW ( j , 1) , 2) , Lcentres (relW( j , 1) , 1) , • go ■ ) ; 

elseif relW (j, 2) <= 0.85 

plot (Lcentres (relW(j , 1) ,2) , Lcentres (relW (j , 1) ,1) , 'y*') ; 
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plot (Lcentres (relW( j , 1) ,2) , Lcentres (relw ( j , 1) ,1) , 'yo») ; 
else 

plot (Lcentres (relw( j , 1) , 2) , Lcentres (relw ( j , 1) , 1) , 1 w* • ) ; 
plot (Lcentres (relW(j , 1) ,2) , Lcentres (relw (j , 1) f l),'wo'); 
end 

hold on; 

end 
end 

end 
end 



function [ivals] =iprof ile (coords, im, pspacing) 

% [ivals] = iprof ile (coords, im) 

% 

% Input (coords) is a matrix where each row contains a 
% 2D coordinate pair (4 values) that specify the end 
% points of a cortical line profile 

% Returns the intensity values at 1 voxel spacing (equiv) 
% If more than one profile is requested, the output rows 
% may be padded with zeros (to the right) . 

% im is the image file and pspacing is the number of points read per 

voxel 

% 

% See also GETCONTOURS, SKELETON, SKELCONTOURS 

X-[] ; 
max 1 enac- 
tor n=l : length (coords ( : , 1) ) , 

x0=coords (n, 1 :2) ; 

xl=coords (n, 3 :4) 

maxlen=max (maxlen, norm (xl-xO) *pspacing) ; 
end 

ivals = zeros ( [length (coords ( : , 1) ) , ceil (maxlen+1) 3 ) ; 

for n=l : length (coords ( : , 1) ) , 
x0=coords (n, 1 : 2) ; 
xl=coords (n, 3:4) ; 
len=norm (xl-xO) *pspacing; 
for r=0 : len, 

x=x0+ (xl-xO) * (r/len) ,- 

ivals (n,r+l) = interp2 (im, x (2) , x (1) , 1 linear ') ; %% swap x &y 
end 
end 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%% 

function [Ledges, Lcentres, num_lam, Id, pcount , W, relw] = 
IPA (coords , ip, std_dev, pi, pspacing, plotnumber) ; 
% analyzes a line profile 
% 

% INPUT: 

% ip is a 2xN matrix containing the interpolated line profile 



WO 02/098292 



37 



PCT/AU02/00731 



% std_dev is the standard deviation of the background noise 

% pi (optional) if pi > 0 then plot the line profile 

% 

% OUTPUT: 

% Id is a list of the intensity drops of the lamina 

% W and relw are a list of the widths of the lamina and the 

coordinates for 

% relative depth of the middle of the laminae 

% num_lam and pcount are the number of lamina and number of peaks 
respectively 

% Ledges and Lcentres are lists of the coordinates for the concavity 
changes and 

% stationary points respectively 
% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%% 

Lcentres= [] ; 
Ledges* [] ; 
c_count=0; 
[x,y] o size (ip) ; 

% find stationary points 

[s_xs, s_ys, s_count, s_ sign, Id, pcount , W, re 1W] = stationary (ip, 
p spacing, std_dev) ; 
if b_xs > 0 

Lcentres=coord (coords, s_xs , pspacing) ; 

end 

count lamina 
num_lam = s_count; 

% find changes in concavity 

[c_xs, c_ys, c_count, c_signj = concavity (ip, std_dev, s_xs, s_sign, 
s_count, pspacing) 

Ledge s=coord (coords, c_xs, pspacing) 
if pi <= 0, return; end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%% 

% plot the line profile 

scalelst_dif f =20; % magnitude of the gradient is small 
scale2nd_dif f = 30; 
figure (plotnumber) 
hold on; 

dip = gradient (ip) ; 

ddip ■ gradient (dip) ; 

plot (l:x, dip.*scalelst_diff , 'r 1 ) ; 

hold on; 

plot (l:x, ddip. *scale2nd__dif f , , g l ) 
hold on; 

plot(l:X,ip(:,l) , «b') ; 

hold on; 
t=l; 

while t<= c_count 

% plot concavities 
if c_count > 0, 
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plot (c_xs (t) , ip(c_xs(t) ,1) , 'r+') ; 
hold on; 

end 
t=t+l; 

end 
t=l; 

while t<= s_count 

% plot stationary points 
if account > 0, 

plot(s_xs(t) , ip(s_xs(t) ,1) , , g+ l ) ; 

end 
t-t+l; 

end 

grid on; 

set the scale for the line profile graphs 
xmin = 0 ; 
xmax = x; 
ymin = -1000; 

ymax = max ( ip ( : , 1 ) ) * 1.4; 
axis ([xmin, xmax, ymin, ymax]); 
xlabel ( 'distance (mm) 1 ) ; 
ylabel ( 1 intensity 1 ) ; 

hold off; 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%% 

function [xs, ys, count, sign, Id,pcount,W, relw] = stationary (ip, 
pspacing, std_ dev) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%% 

if nargin < 2 | nargin > 3 

error ('need a vector containing the intensities along a line 
profile & a search starting point (optional) '); 
end 

W=[] ; 
relW= [] ; 

dip = gradient (ip) ; 
[x,y] = size (ip) ; 
count = 0 ; 
xs(l) = -1; 
ys(l) = -1; 
rsign = 0; 
sign = 0; 
Id=[] ; 

if y < 0, return; end 

[rxs, rys, rcount, rsign] =peaks (ip, std_dev) ; 
temp=l; 

while i<=rcount 

if rcount -i <=2 & rcount -i >= 0 
xs ( count +1) =rxs (i) ; 
ys (count+1) =rys (i) ; 
sign(count+l)=rsign(i) ; 
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count=count +1; 
i-i+1; 

else 

i=i+2 ; 

if rxs (i) -rxs (temp) > 0.5*pspacing 
xs (count +1) srxB (temp) ; 
ys (count+1) =rys (temp) ; 
s ign ( count +l)=rsign( temp) ; 
xs (count+2) =rxs (i~l) ; 
ys (count+2) =rys (i-1) ; 
sign (count +2) =rsign (i-1) ; 
count = count +2 ; 
temp=i; 

else 

if rsign (temp) *ip (rxs (temp) ) <ip (rxs (i) ) 
temp=i ; 

end 

end 

end 

end 

pcount=0 ; 

if xs <0 return; end 

for j =1 : length (xs) 

if sign(j) > 0 

%drop=ip (xs ( j ) ) ; 

if length (xs) ==1 | j==l 

drop= 100-100*ip(l)/ip(xs(j) ) ; 
elseif j<length(xs) 

drop=100-100*abs ( (ip (xs ( j - 
l))+ip(xs(j+l)))/2)/ip(xs(j) ) ; 
else 

drop=100-100*ip(xs(j-l))/ip(xs(j) ) ; 

end 

Id= [Id; j ,drop] ; 
pcount«pcount+l ; 

end 

end 

for k=l : length (xs ) 

if sign(k) > 0 

if length (xs) ==1 

width=0 . 5 * length (ip) /length (ip) ; 
elseif k==l 

width=0 . 5* (xs (k+1) -1) /length (ip) ; 
elseif k<length(xs) 

width=0 . 5* (xs (k+1) -xs (k-1) ) /length (ip) ; 

else 

width=(0.5*(xs(k) -xs (k-1) ) + length (ip)- 

xs (k) ) /length (ip) ; 

end 

W= [W;k, width] ; 

end 

end 

for m=l : length (xs) 

if sign(m) > 0 

rw=xs (m) /length (ip) ; 
relw= [relW;m # rw] ; 

end 

end 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%% 

function [xs, ys, count, sign] = concavity (ip, std_dev, pspacing);% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%% 

if nargin < 3 | nargin > 4 

error <* need a vector containing the intensities along a line 
profile & a search starting point (optional)'); 
end 

dip = gradient (ip) ; 

count = 0 ; 

xs (1) = -1; 

ys(l) = -1; 

rsign = 0 ; 

sign = 0; 

[rxs, rys, rcount, rsign] =peaks (dip, std_dev) ; 

temp=l; 

i = l; 

while i<=rcount 

if rcount -i <=2 & rcount -i >= 0 
xs (count +1) =rxs (i) ; 
ys (count +1) =rys (i) ; 
sign (count+1) =rsign(i) ; 
count =count +1; 
i=i+l; 

else 

i=i+2 ; 

if rxs (i) -rxs (temp) > 0.5*pspacing 
xs ( count +1) =rxs (temp) ; 
ys ( count +1) =rys (temp) ; 
sign (count +1) =rsign (temp) ; 
xs (count +2) =rxs (i-1) ; 
ys (count +2) =rys (i-1) ; 
sign (count+2) =rsign (i-1) ; 
count e count +2 ; 
temp=i ; 

else 

if rsign (temp) *dip (rxs (temp) ) <dip (rxs (i) ) 
temp=i ; 

end 

end 

end 

end 

p count =0; 

if xs <0 return; end 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%% 

function [rxs, rys, rcount, rsign] = peaks (ip, std_dev);% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%% 

if nargin < 1 



i 
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error ('need a vector containing the intensities along a line 
profile ' ) ; 
end 

dip = gradient (ip) ; 
[finish, y] = size(ip); 
r count = 0 ; 
i = 1; 

temp = i ; 
rsign = 0; 

while i<finish 

if (dip(i)*dip(i+l) <= 0) 
temp = i; 

while ( (dip(i) *dip(i+l) ) ==0) 
i = i + 1; 
if i==finish 
return; 

end 

end 

if (i mm temp) & i < finish 
i - i + 1; 

end 

if ( (dip (temp) *dip(i) ) < 0 ) | (i « finish & temp < 

finish) 

if i-temp ==1 

if dip (temp) < 0 
ppos=temp; 

else 

ppos=i ; 

end 

else 

ppos= round ( (i+temp) /2) ; 

end 

rxs (r count +1) = ppos; 
rys (rcount+1) = ip(ppos); 
rcount = rcount + 1; 
if rcount > 0 

if (dip (temp) > 0) & rcount > 0 

rsign (rcount) = 1; 
elseif dip (temp) < 0 

rsign (rcount ) = -1; 

end 

end 

end 

else 

i = i+1 ; 

end 

end 

if rcount > 0 

if rsign (rcount ) == -1 

if ip (length (ip) ) > ip (rxs (rcount ) ) 
rxs (rcount+1) =length(ip) ; 
rys (rcount+1) =ip (length (ip) ) ; 
rcount=rcount+l ; 
rsign (rcount) =1; 

end 

end 

end 

if rcount < 1; 

rxs(l)=-l; 
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rys (1) =-1; 

end 



function [Points] =coord (coords , xs, pspacing) 

Points= [] ; 
xO=coords (1:2) ; 
Xl=COOrds (3:4) ; 
len=norm(xl~xO) *pspacing; 



for r=l : length (xs) , 

x= (xs (r) /len) * (xl-xO) + xO; 
%f igure(l) ; 

%plot (x(l,2) , x(l,l) , 'g-o') ; 
%hold on; 

Points = [Points ; x] / 

end 



Computer code used for 3-D analysis 

function V2 = anzload(anzfile,x,y,z,anzprecision,machineformat) 
if nargin = 5 

machineformat = T; 
elseif nargin > 6 | nargin < 5 

si = sprintfCneed 5 (or an optional 6th) arguments:\n'); 

s2 = sprintf('l: the filename of an Analyze image (.img)\n'); 

s3 = sprintf('2,3,4: voxel dimensions of the imageW); 

s4 = sprintfCS: datatype used in the image\n f ); 

s5 = sprintf( ? 6: (optional) format of the data\n'); 

s = Strcat(sl,s2,s3,s4,s5); 

error(s); 

end 

fid = fopen(anzfile, 'r\ machineformat); 
iffid = -l 

s = sprintf('cannot open file \"%s\"\n f , anzfile); 
error(s); 

end 

[VOL,count] = firead(fid, [x*y*z], anzprecision); 
fclose(fid); 

V2 = reshape(VOL, [x,y,z]); 



function [start, stop] = boundary j>ts(lp,c_xs, c_ys, c_count,spacing, grey_range, 
c_sign, std_dev, VOX); 

if nargin < 1, error('need a line profile vector 1 ); end 

dip = gradient(lp); 

[x,y]=size(lp); 

magnitudes = zeros(c_count,l); 
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start = -1; 

stop = -1; 

max__peak = -1; 

maxmag = [0,0]; 

if c__count <= 1, return; end 

for j=l:c_count 

if c_ys(j) > max_peak 
max_peak = c_ys(j); 

end 

end 

if max_peak <= 0, return; end 
for i=l :c_count 

el = c_ys(i) / max_peak ; 

e2 = (x - c_xs(i)) / x; 

tmp = abs(grey_range(l) - lp(c_xs(i))); 

if lp(c_xs(i)) > grey_range(l) 

e3 = (1 - tmp/(grey_range(2)-grey_range(l))); 

else 

e3 = (1 - 4*tmp/(grey_range(2)-grey_range(l))); 

end 

ife3<0 
e3 = 0; 

end 

if c_sign(i) > 0 

magnitudes(i) = sqrt(el A 2 + e2 A 2 + e3 A 2); 

else 

magnitudes(i) = 0; 

end 

end 

for i=l:c_count-l 
if i < c_count 

if magnitudes(i)> max_mag(2) 
maxmag = [i,magnitudes(i)]; 

end 

end 

end 

base = x; 
for i=l:x 

iflp(i)<2 
base = i; 
break; 

end 

end 

if max_mag(l) <= c_count - 2 

e = max_mag(l)+ 2; 
elseif max_mag(l) >c_count - 2 

e = maxmag(l); 

end 

if max mag(l) <= 2 
b = max_mag(l); 
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elseif maxmag(l) > 2 
b = max_mag(l) - 2; 

end 

if b <= 0 

start = -1; 
stop = -1; 
return; 

end 

temp = c_ys(b); 
start=c_xs(b); 
while b < e 

if c_sign(b) > 0 & c_ys(b) > temp & abs(c_xs(b)-c_xs(max_mag(l))) < 
max(VOX)/(2*spacing) 
temp = c_ys(b); 
start = c_xs(b); 

end 

b = b + 1; 

end 

if base > start 
start = -1; 
break; 

end 

magnitudes = zeros(c_count,l); 
for i=l:c_count 

el = 0; 

e2 = 0; 

e3 = 0; 

if ((c_xs(i)-start)*spacing < 5.5) & ((c_xs(i)-start)*spacing > 2) 
el = c_ys(i) / max_peak; 

e2 = l-abs((start+(4/spacing) - c_xs(i)) / (start+(4/spacing))); 
if lp(c_xs(i)) < grey_range(2) 

e3 = (l-abs(grey_range(2)-lp(c_xs(i)))/(grey_range(2)-grey_range(l))); 

else 

e3 = (l-4*abs(grey - range(2)-lp(c_xs(i)))/(grey_range(2)- 
grey_range(l))); 
end 

end 

if el <= 0 | e2 <= 0 | e3 <= 0 
el=0; 
e2 = 0; 
e3=0; 

end 

if c_sign(i)>0 

if c_xs(i)~= start 

magnitudes(i) = sqrt(el A 2 + e2 A 2 + e3 A 2); 

else 

magnitudes(i) = 0; 

end 

else 

magnitudes(i) = 0; 
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end 

end 

maxmag = [0,0]; 
for i=l:c_count 

if magnitudes(i)> max_mag(2) 

max_mag = [i,magnitudes(i)]; 

end 

end 

if maxmag(l) — 0 

stop = c_xs(max_mag(l )); 

else 

stop = -1; 

end 

if base > start, stop = -1; 

break; 

end 



function [R,G,B] = colour_picker(drops, max drop, widths, numlam, type); 
if nargin — 5 

error('need drop, max drop, width, number of lamina and type of colouring'); 

end 

if existCDEFINE^SCHEME 1 ) — 1 
colour_scheme 

end 

switch type 

case { NUM_LAM } 
if num_lam <= 0 

R = 255; 

G = 255; 

B = 255; 
elseif num_lam = 1 

R = 0; 

G = 0; 

B = 255; 
elseif num lam = 2 

R = 0; 

G = 255; 

B = 0; 

else 

R = 255; 
G = 0; 
B = 0; 

end 

case { FIRST_DROP } 
ifdrops(l)>0 

h = 0;s= 1; v = 255; 
s = drops(l)/max_drop; 
[R,G,B] = hsv2rgb([h,s,v]); 
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round(R); round(G); round(B); 
elseif drops(l) < 0 
R = 0; 
G = 0; 
B = 255; 

else 

R = 255; 
G = 255; 
B = 255; 

end 
otherwise 

if num_lam > 1 

maxd = max(drops); 
for i=l:num_lam 

if drops(i) = maxd, break; end 

end 

h = 0; s= 1; v = 255; 

s = drops(i) / 1 0; 

[R,G,B] = hsv2rgb([h,s,v]); 

R « round(R); G = round(G); B = round(B); 
elseif num_lam = 1 

h = 0;s = l; v = 255; 

s = drops(l)/ 10; 

[R,G,B] = hsv2rgb([h,s,v]); 

R = round(R); G = round(G); B = round(B); 
elseif numlam = 0 

R = 255;G = 255;B = 255; 

else 

R = 0;G = 0;B = 255; 

end 

end 



global DEFINE SCHEME; 

global NUM__LAM FIRST_DROP MAX^DROP FIRST WIDTH MAX_WIDTH; 

NUM_LAM = 1; 

FIRST_DROP =2; 

MAX_DROP = 3; 

FIRSTWIDTH = 4; 

MAXJWTDTH =5; 



function [xs, ys, count, sign] = concavity(lp, std_dev); 
if nargin < 2 | nargin > 3 

error( f need a vector containing the intensities along a line profile & a search 
starting point (optional) 1 ); 
end 

dip = gradient(lp); 
[x,y] = size(lp); 

[xs, ys, count, sign] = peaks(dlp, 1, y, std_dev); 
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function [drops, widths, numjam] = lpa(lp, spacing, num, std dev, pi, tissuetype, 
grey_range, VOX); 
if nargin ~= 8 

error('wrong number of input arguments 1 ); 

end 

OK=l; 
BAD = 0; 
drops = 0; 
widths = [0,0,0]; 
tag = OK; 
app_s = "; 

if nargin = 2, pi = 0; end 
[x,y] = size(lp); 
xmin = 0; 
xmax = lp(x,l); 
ymin = -5; 

ymax = grey_range(2)*2; 

[c_xs, c_ys, c_count, c__sign] = concavity(Ip(:,2)', std_dev*spacing A 4); 
[start,stop] = boundary_pts(lp(:,2),c_xs, c_ys, c_count, spacing, grey_range,c_sign, 
std_dev*spacing, VOX); 
if start <0 | stop<0 
tag = BAD; 

app_s = strcat(app_s,' [cannot find boundaries]*); 
numlam = -1; 

end 

[s__xs, s_ys, s_count, s_sign] = stationary(lp(:,2) , , start, stop, tissue_type,std_dev, 
spacing); 

num_lam = round(s_count/2); 

if ((s_count = 0) | (s_xs(s_count) < start)) & (tag ~= BAD) 
app_s = strcat(app_s,' [no stationary pts] 1 ); 

end 

ind = 1 ; 
iftag = BAD 

drops(l) = -1; 
num_lam = -1; 
elseif s_count > 1 & start >= 0 
for j=l :(s_count- 1 ) 

tmp = s_ys(j) - s_ys(j+l); 
iftmp>0 

if std_dev — 0 

drops(ind) = tmp / std_dev; 

else 

drops(ind) = tmp; 

end 

for k=l :(c_count-l) 

if c_xs(k) < s_xs(j) & c_xs(k+l) > s_xs(j) 
ca = c__xs(k)-s_xs(j); 
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cb = (c_xs(k)-s_xs(j))/(s_xs(j+l) - s_xs(j)); 
cc = (s__xs(j+l) - s_xs(j)); 
widths(ind,:) = [ca,cb,cc]; 

else 

ca = (s_xs(j+l) - s_xs(j))/2; 
cb = 0; 

cc = (s_xs(j+l) - s_xs(j)); 
widths(ind,:) — [ca,cb,0]; 

end 

widths(ind,:) = widths(ind,:).*spacing; 

end 

ind = ind + 1 ; 

end 

end 

end 

if pi <= 0, return; end 
scale lst_diff = 20; 
scale2nd_diff = 100; 
if tag ~= BAD 

dip = gradient(lp(:,2)); 

ddlp = gradient(dlp); 

plot(lp(:,l), dlp.*scalelst_diff,V); 

hold on; 

plot(lp(:, 1 ), ddlp.*scale2nd_diflf, , g*); 

end 

plot(lp(: 5 l),lp(:,2),V); 
hold on; 

if s_count > 0, plot(lp(s_xs(l),l), s^OVr*'); end 
if s_count > 1, plot(lp(s_xs(2),l), s_ys(2),'g+'); end 
if start >0& stop >0 

plot(start*spacing, lp(start,2), Ic**); 

plot(stop*spacing, lp(stop,2), Tc* 1 ); 

end 

grid on; 

axis([xmin, xmax, ymin, ymax]); 
iftag = OK 

foo = strcat(*line profile (max lamination drop = num2str(drops(l)), y ,app_s); 
title(foo); 

else 

foo = streamline profile (DUD)\app_s); 
title(foo); 

end 

xlabel('distance (mm) 1 ); 
ylabelCintensity'); 
hold off; 
pause; 



function [newx, newy, newz] = normalize(oldx, oldy, oldz); 

if nargin — 3, errorCnormalize works on 3D points/vectors only 1 ); end 
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ssum = sqrt(abs(oldx) A 2 + abs(oldy) A 2 + abs(oldz) A 2); 
if ssum = 0, 

waming( f cannot find normal, possible degenerate triangle.'); 

newx = 0; 
newy = 0; 
newz = 0; 

else 

newx = oldx / ssum; 
newy = oldy / ssum; 
newz = oldz / ssum; 

end 



function [rxs, rys, rcount, rsign] = peaks(lp, start, finish, std_dev); 
if nargin < 1 

error( f need a vector containing the intensities along a line profile'); 

end 

dip = gradient(lp); 
[x,y] = size(lp); 
count = 0; 
rcount = 0; 
i = start; 
temp = i; 
sign = 0; 
rsign = 0; 
while i<finish 

if(dlp(i)*dlp(i+l)<=0) 
temp = i; 

while (abs(dlp(i)*dlp(i+l)) <= std_dev)& (i<finish-l) 
i = i+ 1; 

end 

if i = temp & i < finish 
i = i + 1; 

end 

if ((dlp(temp)*dlp(i)) < 0 & dlp(temp) ~= dlp(i)) | i == finish- 1 
rxs(rcount+l) = round((i + temp)/2); 
rys(rcount+l) = lp(rxs(rcount+l)); 
rcoimt = rcount + 1 ; 
if rcount > 0 

if dlp(temp) > 0 & rcount > 0 

rsign(rcount) = 1 ; 
elseif dlp(temp) < 0 
rsign(rcount) = -1; 

end 

end 

end 



else 

i = i+l; 
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end 

end 

if rcount <= 1; 
rxs(l)=-l; 
rys(l)=-l; 

end 



% quickrun.m 

clear all; 
close; 

format long 
tic 

if exist( , DEFINE_SCHEME , ) — 1 
colour_scheme 

end 

ot = MAX_DROP; 
nl = sprintfCW); 

dispCPlease enter the name of a mesh file. 1 ); 

fin = deblank(input('mesh file: \W)); 

[path, base, ext, ver] = fileparts(fin); 

s = sprintfCoutput will be in file \ M %s.off\"\n\ base); 

disp(s); 

fout = strcat(base,\out'); 
fid = fopen(fout,'w'); 

disp('Depth (mm) the line profile should reach into the cortex? 1 ); 

depth = input('depth: '); 

disp(nl); 

dispCNumber of sample points to interpolate? 1 ); 

samples = input( f samples: '); 

disp(nl); 

disp( , Standard deviation of background noise? 1 ); 

disp('(Enter the value 1, to ignore signifigance rel. to background noise.)'); 

std_dev = input('std. deviation: '); 

disp(nl); 

dispCPlease enter the lower limit of the grey matter intensity range? 1 ); 

low = input('intensity value: '); 

disp(nl); 

dispCPlease enter the upper limit of the grey matter intensity range?*); 
high = input('intensity value: '); 
grey_range = [low,high]; 
disp(nl); 

disp('What type of tissue analysis do you wish to perform?'); 
disp('choices: '); 

disp( f,, g" : grey matter analysis 1 ); 
dispC'w" : white matter analysis'); 
tissue type = input('tissue type: ','s'); 
disp(nl); 

dispCWhat kind of analysis do you want to perform?'); 
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disp('choices: *); 

s=sprintf( m %d" : colour coded for number of laminations^NUM^LAM); 
disp(s); 

^sprintfrTod" : colour coded for first intensity drop , ,FIRST_DROP); 
disp(s); 

disp( ,f, 0" : colour coded for maximum intensity drop 1 ); 
ot = input('analysis type: '); 
hfile = strcat(path, V, base, '.img 1 ); 
switch lower(ext) 
case {'.wfr 1 } 

[num_tri angles, num_points, triangles, points] = read_emse(fin,hfile); 
otherwise 

s = sprintf^unknown mesh type \ ,, %sV'\n , ,ext); 
error(s); 

end 

[DIM VOX SCALE TYPE OFFSET ORIGIN DESCRIP] = spm_hread(hfile); 
fprint^fid,^ %d\n*, num_triangles); 
pi = zeros(l,3); 
p2 = zeros(l,3); 
p3 =zeros(l,3); 

j-i; 

for i=l : 1 :num_tri angles 

pi = points(triangles(i,l)+l,:); 
p2 = points(triangles(i,2)+l,:); 
p3 = points(triangles(i,3)+l,:); 

pl(l) = pl(l)*VOX(l); 
P 2(l)«p2(l)*VOX(l); 
P 3(l)=p3(l)*VOX(l); 
pl(2) = pl(2)*VOX(2); 
P 2(2)=p2(2)*VOX(2); 
P 3(2)=p3(2)*VOX(2); 
pl(3) = pl(3)*VOX(3); 
P 2(3)=p2(3)*VOX(3); 
P 3(3)=p3(3)*VOX(3); 

[cx, cy, cz] = tricenter(pl, p2, p3); 
[nx, ny, nz] = trinormal(pl, p2, p3); 
[nx, ny, nz] = normalize(nx,ny,nz); 
cx = cx / VOX(l); nx = nx / VOX(l); 
cy = cy / VOX(2); ny = ny / VOX(2); 
cz = cz / VOX(3); nz = nz / VOX(3); 
if[nx,ny,nz] = [0,0,0]; 

dispftriangle number:*); 

disp(i); 

else 

fprintfitfid/s %f %f %f d %f %f %f\n , ,cx,cy,cz, nx,ny,nz); 

j=j + i; 

end 

end 
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fclose(fid); 

snc(base, depth, samples, stdjiev, ot, 0, tissue_type, grey_range); 
toe 



function [num_tri, num_jjoints, triangles, points] = read_emse(fhame, hname); 
if nargin — 2 

error( f need meshfile and header file 1 ); 

end 

fid = fopen(fhame, 'it 1 ); 
iffid = -l 

s = sprintfCcannot open file \"%s\"\n\ fiiame); 
error(s); 

end 

[path, base, ext, ver] = fileparts(fhame); 

[DIM VOX SCALE TYPE OFFSET ORIGIN DESCRIP] = spm Jiread(hname); 

output_file = strcat(base, '.off); 

tri_file = strcat(base, f .tri'); 

fout — fopenfautput^file/w'); 

for i=l:l:3, fgetl(fid); end 

str = sprintfCgrep -c v %s\ fhame); 

[s,w] = unix(str); 

num_points = str2num(w) 

points = zeros(num_points, 3); 

str = sprintfTgrep t %s | sed Y's/t/S/gV > %s\ fiiame, tri^file); 
unix(str); 

str = sprintflCgrep -c t %s', fiiame); 
[s,w] = unix(str); 
num_tri = str2num(w) 
triangles = zeros(num_tri, 3); 

fprintf(fout,'OIT\n 0 /od %d %d\n', num_points, num_tri, 0); 

i-i; 

p_ind=l; t_ind=l; 
line = fgetl(fid); 
while line — -1, 

P = sscanf(line, V %f %f %f)'; 
if ~isempty(P) 

a = P(l);b = P(2);c = P(3); 
points(p__ind, :)=[a,b,c] ; 
p_ind = p_ind + 1; 

fprintf(fout,«°/of %f %f\n', a*VOX(l), b*VOX(2), c*VOX(3)); 

else 

T = sscanf(line, »t %i %i %i*)'; 
if~isempty(T) 

triangles(t Jnd, :)=T; 

t_ind = t_ind + 1 ; 

end 

end 
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line = fgetl(fid); 

end 

fclose(fid); 



function snc(base, depth, N, std_dev, ot, pi, tissue_type, grey_range); 
if nargin < 7 | nargin > 8; 

errorCwrong number of arguments 1 ); 

end 

if existCDEFrNE.SCHEME 1 ) — 1 
colour_scheme 

end 

output_file = strcat(base, '.off); 
intfile = strcat(base, '.out'); 
tri file = strcat(base, '.tri'); 
img file = strcat(base, , .img'); 

[DIM VOX SCALE TYPE OFFSET ORIGIN DESCRIP] = spmjiread(img_file); 

vol_data = anzload(img_file, DIM(l), DIM(2), DIM(3), 'uint8\ T); 

fout = fopenCoutpu^file/a'); 

fid = fopen(int_file, V); 

ft = fopenttrLfile/r'); 

line = fgetl(fid); 

num_tri angles = sscanf(line, f # %d*); 

k = depth/N; 

for i=l :num_tri angles 

line = fgetl(fid); 

tri = fgetl(ft); 

A = sscanfUine/s %f %f %f d %f %f %f); 

xO = A(l); yO = A(2); zO = A(3); 

nx = A(4); ny = A(5); nz = A(6); 

dx = k*nx; dy = k*ny; dz = k*nz; 

x = (xO + (-N/2:N-l)*dx); 

y = (y0 + (-N/2:N-l)*dy); 

z = (zO + (-N/2:N-l)*dz); 

scalev = (l:3/2*N)*k; 

lp = inteip3(vol_data,x,y,z, , *Iinear'); 

newlp = [ scalev',lp']; 

[drops, widths, num_lam] = lpa(newlp, k, i, std_dev, pi, tissue__type, grey_jange, 
VOX); 

max drop = grey_range(2) - grey_range(l); 

[R, G, B] = colour_picker(drops, max drop, widths, num_lam, ot); 
fprintf(fout, ,0 /os %d %d %d\n\ tri, R, G, B); 

end 

fclose(fid); 
delete(tri_file); 



function [xs, ys, count, sign] = stationary(lp, agbound, gwbound, tissue_type, std_dev, 
spacing); 
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if nargin < 5 | nargin > 6 

error( f need a vector containing the intensities along a line profile & a search 
starting point (optional)'); 
end 

[x,y] = size(lp); 
count = 0; 
xs(l) = -l; 
ys(l) = -l; 
rsign = 0; 
sign = 0; 

if tissue_type = f g' 

start = agbound; 

stop = gwbound; 
elseif tissue_Jype — W 

start = gwbound; 

stop = y; 

else 

errorCunknown tissue type 1 ); 

end 

if stop < 0, return; end 

[rxs, rys, rcount, rsign]=peaks(lp, start, stop, std_dev*spacing A 2); 
if rcount > 2 

if rys(rcount) > rys(rcount -1) 
rcount=rcount - 1 ; 

end 

end 

for i=l:2:rcount-l 

if lp(rxs(i))-lp(rxs(i+l)) > 2 * std_dev 
xs(count+l)=rxs(i); 
xs(count+2)=rxs(i+ 1 ) ; 
ys(count+ 1 )=rys(i); 
ys(count+2)=rys(i+ 1 ) ; 
count=count+2; 
sign(count+ 1 )=rsign(i) ; 
si gn(count+2)=rsign(i+ 1 ) ; 

end 

end 



% testrun.m 

clear all; 
close; 

format long 
tic 

if existCDEFINE^SCHEMES') — 1 
colour_scheme 

end 

ot = MAX_DROP; 
nl = sprintf^); 
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disp(Tlease enter the name of a mesh file. 1 ); 

fin = which(deblank(input( , mesh file: 7s 1 ))); 

[path, base, ext, ver] = fileparts(fin); 

s = sprintfCoutput will be in file V , %s.offV , \n , > base); 

disp(s); 

fout = strcatCbase/.out 1 ); 
fid = fopen(fout,Sv f ); 

disp(T>epth (mm) the line profile should reach into the cortex?'); 

depth = input( f depth: '); 

disp(nl); 

dispCNumber of sample points to interpolate? 1 ); 

samples = input('samples: f ); 

disp(nl); 

disp('Standard deviation of background noise? 1 ); 

disp('(Enter the value 1, to ignore signifigance rel. to background noise.) 1 ); 

std_dev = input('std. deviation: '); 

disp(nl); 

disp('Please enter the lower limit of the grey matter intensity range?'); 

low = input('intensity value: '); 

disp(nl); 

disp(Tlease enter the upper limit of the grey matter intensity range?'); 

high = input('intensity value: '); 

disp(nl); 

dispCWhat type of tissue analysis do you wish to perform?'); 
disp('choices: f ); 

disp( m g" : grey matter analysis'); 
disp("'w" : white matter analysis'); 
tissue_type = input('tissue type: ','s'); 
disp(nl); 

grey_range = [low,high]; 
hfile = strcat(path, V, base, '.img'); 
switch lower(ext) 
case {'.wfr'} 

[nuinjri angles, num__points, triangles, points] = read_emse(fin,hfile); 
otherwise 

s = sprintf('unknown mesh type Y'%s\"\n',ext); 
error(s); 

end 

[DIM VOX SCALE TYPE OFFSET ORIGIN DESCRIP] = spm_hread(hfile); 
fprintf(fid,'# %d\n\ num_triangles); 
pi =zeros(l,3); 
p2 = zeros(l,3); 
p3 = zeros(l,3); 

j=i; 

for i=l : 1 :num_tri angles 

pi =points(triangles(i,l)+l,:); 
p2 =points(triangles(i,2)+l,:); 
p3 =points(triangles(i,3)+l,:); 
pl(l) = pl(l)*VOX(l); 
P 2(l) = p2(l)*VOX(l); 
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P 3(l) = p3(l)*VOX(l); 

pl(2) = pl(2)*VOX(2); 

P 2(2) = p2(2)* VOX(2); 

p3(2) = p3(2)* VOX(2); 

pl(3) = pl(3)*VOX(3); 

P 2(3) = p2(3)* VOX(3); 

p3(3) = p3(3)* VOX(3); 

[cx, cy, cz] = tricenter(pl, p2, p3); 

[nx, ny, nz] = trinormal(pl, p2, p3); 

[nx, ny, nz] = normalize(nx,ny,nz); 

cx = cx / VOX(l); nx = nx / VOX(l); 

cy = cy / VOX(2); ny = ny / VOX(2); 

cz = cz / VOX(3); nz = nz / VOX(3); 

if [nx,ny,nz] = [0,0,0]; 

disp('triangle number: 1 ); 

disp(i); 

else 

rprintf(fld,'s %f %f %f d %f %f %f\n',cx,cy,cz, nx,ny,nz); 

j-j + i; 

end 

end 

fclose(fid); 

snc(base, depth, samples, std_dev, ot, 1, tissue_type, grey_range); 
toe 



function [x,y,z] = tricenter(pl, p2, p3); 
x = (pl(l)+p2(l)+p3(l))/3; 
y = (pl(2)+p2(2)+p3(2))/3; 
z = (pl(3)+p2(3)+p3(3))/3; 



function [cx,cy,czj = trinormal(pl, p2, p3); 

if nargin ~= 3, error( t need three 3D points as input'); end 

ax = pl(l)- P 2(l); 

ay = pl(2)-p2(2); 

az = pl(3)-p2(3); 

bx = p3(l)-p2(l); 

by = p3(2)-p2(2); 

bz = P 3(3)-p2(3); 

cx '— ay*bz - az*by; 

cy = az*bx - ax*bz; 

cz = ax*by - ay*bx; 
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THE CLAIMS DEFINING THE INVENTION ARE AS FOLLOWS: 

1. A method of mapping a property of a three dimensional 
object, said method comprising the steps of 
mapping the property for at least a portion of the 
object , 

providing a line or region which defines intersections 
with part of said mapped portion, and 
displaying the property for the intersections. 

2. A method of mapping a property of a three dimensional 
object, said method comprising the steps of 
mapping the property for a plurality of slices within 
the object, 

providing a line or region which defines intersections 
with part of at least some of the plurality of slices, 
and 

displaying the property for the intersections. 

3. A method as defined in claims 1 or 2 wherein the 
property relates to nuclear magnetic resonances. 

4. A method as defined in claim 2 wherein mapping of the 
property for the plurality of slices results in a 
plurality of nuclear magnetic resonance images. 

5. A method as defined in any one of the preceding claims 
wherein the line or region is a line. 



6. A method as defined in claim 5 wherein the line is a 
plurality of lines. 
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7. A method as defined in claim 6 wherein each of the 

lines is substantially perpendicular to the surface of 
the three-dimensional obj ect . 



8. A method as defined in any one of the preceding claims 
which further comprises the step of approximating the 
shape of the object by a three dimensional lattice of 
two-dimensional forms created by a first computer 
routine . 



9. A method as defined in claim 8 wherein the lattice is 
a mesh. 



10. A method as defined in any one of claims 8 to 9 
wherein each of the lines is a line substantially 
perpendicular to one of the two-dimensional forms. 

11. A method as defined in claim 2 wherein the slices are 
selected using a second computer routine. 

12. A method as defined in any one of claims 3 to 11 
wherein a third computer routine maps nuclear magnetic 
resonance values for the intersections by displaying 
an array of values . 

13. A method as defined in claim 12 which further 
comprises the step of analysing the array of values 
using a fourth computer routine. 



14. A method as defined in claim 13 wherein the fourth 
computer routine analyses numerical properties. 
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15. A method as defined in any one of claims 12 to 14 
wherein the step of analysing the array of values 
comprises differentiation. 

16. A method as defined in claim 15 wherein analysing the 
array of values comprises determining the positions of 
the maxima, minima and/or zero-points of the array of 
values and the derivatives of the array of values. 



17. A map produced by a method of mapping a property of a 
three dimensional object, said method comprising the 
steps of 

mapping the property for at least a portion of the 
object, 

providing a line or region which defines intersections 
with part of said mapped portion, and 
displaying the property for the intersections. 

18. A map produced by a method of mapping a property of 
a three dimensional object, said method comprising the 
steps of 

mapping the property for a plurality of slices within 
the object, 

providing a line or region which defines intersections 
with part of at least some of the plurality of slices, 
and 

displaying the property for the intersections. 

19. A method of mapping a property of an object, said 
method comprising the steps of 

mapping the property for a slice within the object, 
providing a line within the slice, and 
displaying the property for the line. 
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20. A method as defined in claim 19 wherein the property 
wherein relates to nuclear magnetic resonances. 

21. A method as defined in any one of claims 19 to 20 
wherein mapping of the property for the slice results 
in a nuclear magnetic resonance image. 

22. A method as defined in any one of claims 19 to 21 , 
wherein the line is a plurality of lines. 

23. A method as defined in any one of claims 19 to 22 
wherein a fifth computer routine maps nuclear magnetic 
resonance values for the intersections by displaying 
an array of values. 

24. A method as defined in claim 23 which further 
comprises the step of analysing the array of values 
using a sixth computer routine. 

25. A method as defined in claim 24 wherein the fifth 
computer routine analyses numerical properties. 

26. A method as defined in any one of claims 23 to 25 
wherein the step of analysing the array of values 
comprises differentiation. 

27. A method as defined in claim 26 wherein the step of 
analysing the array of values comprises determining 
the positions of the maxima, minima and/or zero-points 
of the array of values and the derivatives of the 
array of values . 
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28. A map produced by a method of mapping a property of 
an object, said method comprising the steps of 
mapping the property for a slice within the object , 
providing a line within the slice, and 
displaying the property for the line. 
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